clc;clear;

syms z x y l1 l2 l3 theta3 theta4
f1 = (z-l2*cos(theta3)-l3*cos(theta4))^2
f2 = (sqrt(x^2+y^2)+l2*sin(theta3)-l3*sin(theta4))^2
f = f1+f2
fe = expand(f1+f2)
simplify(fe)
S = solve(f==l1^2,theta3, 'Real', true)